################################################################################
# 04-figure02.R
# Paper:    Bottom Up? Top Down? Determinants of Issue-Attention in State 
#           Politics.
# Journal:  Journal of Politics.
# Authors:  Andreu Casas, Oscar Stuhler, Julia Payson, Joshua A. Tucker, 
#           Richard Bonneau, and Jonathan Nagler.
# Purpose:  This script fits the main models for the results presented in
#           sections 5.1 and 5.2 of the paper. It then generates Figure 2.
# Note:     Models were fit on NYU's HPC cluster. Instances used 40 CPU cores
#           and 100 GB of memory. With this setup, fitting both models took
#           half a day. In this script, we comment out the portion where we fit 
#           the models and instead load the fitted model objects from disk.
# Data In:  
#           Two datasets in which we have group-day-issue-level information about
#           how much attention each group of analysis paid to each day to each 
#           topic we study:
#           /data/group-day-issue-level-dataset-01-2018
#           /data/group-day-issue-level-dataset-01-2021
# Data Out:
#           1. Two estimated VAR models, one for 2018 and one for 2021.
#              /models/MODEL_I_2018.Rdata
#              /models/MODEL_I_2021.Rdata
#           2. Two IRF objects with the impulse response estimates based on the
#              fitted models.
#              /models/MODEL_I_IRFs_2018.Rdata
#              /models/MODEL_I_IRFs_2021.Rdata
#              
#           3. Figure 2 of the main paper:
#              - ./figures/figure02.png
################################################################################


#===============================================================================
# PACKAGES
#===============================================================================
library(dplyr)
library(tidyr)
library(vars)
library(ggplot2)
library(ggh4x)
library(boot)


#===============================================================================
# FITTING MODELS
#===============================================================================
#
####### Defines years variable
#years <- c("2021", "2018")
#
####### For each year, run the models with the same parameters
#for(year in years){
#  
#  cat(year, "\n")
#  #-----------------------------------------------------------------------------
#  # DATA
#  #-----------------------------------------------------------------------------
#  # Model data
#  load_file <- paste0("./data/group-day-issue-level-dataset-01-", year, ".csv")
#  maindb <- read.csv(load_file)
#  
#  # Sort:
#  maindb <- maindb[order(maindb$IssueState, maindb$Date),]
#  
#  #-----------------------------------------------------------------------------
#  # MAIN
#  #-----------------------------------------------------------------------------
#  # Create a formula object that will facilitate transforming the topic variable 
#  # from factor to dummies.
#  if(year == "2021"){
#    variables <- c("national_legislators", "national_media", "state_legislators", 
#                   "state_media", "state_partisans", "state_random_partisans", "Trump", "IssueState")
#  }
#  if(year == "2018"){
#    variables <- c("national_legislators", "national_media", "state_legislators", 
#                   "state_media", "state_partisans", "Trump", "IssueState")
#  }
#  mformula <- formula(paste0("~", paste0(variables, collapse = " + ")))
#  model_data <- model.matrix(mformula, maindb[, variables])
#  model_data <- model_data[, 2:ncol(model_data)] # remove intercept
#  
#  # Splitting the covariates of interest from the issue dummy variables
#  X_endogenous <- model_data[, variables[which(!variables %in% c("IssueState"))]]
#  X_exogenous <- model_data[
#    , which(!colnames(model_data) %in% 
#              variables[which(!variables %in% c("IssueState"))])]
#  
#  # Define lags, number of bootstraps, and set seed
#  p <- 5
#  runs <- 300
#  seed <- 78223
#  set.seed(seed)
#  
#  # Run the VAR model and estimate IRFs
#  var_model_merged <- VAR(y = X_endogenous, p = p, exogen = X_exogenous)
#  var_irfs_cum_merged <- irf(var_model_merged, n.ahead = 15, ortho = F,
#                             cumulative = TRUE, runs = runs)
#  
#  # Save the estimated VAR model and IRFs.
#  save(var_model_merged, file = paste0("./models/MODEL_I_", year, ".Rdata"))
#  save(var_irfs_cum_merged, file = paste0("./models/MODEL_I_IRFs_", year, ".Rdata"))
#}



#===============================================================================
# VISUALIZE THE MODELS IN FIGURE 2 USING THE IRFs
#===============================================================================

# Load the figure code
source(paste0("./code/00-functions.R"))

# Define which type of effect to show
which_effects <- "structural"

# Load the IRFs for Model I
year <- "2021"
model_name <- paste0("./models/MODEL_I_IRFs_", year, ".Rdata")
print(load(model_name))
var_irfs_cum_merged_2021 <- var_irfs_cum_merged
year <- "2018"
model_name <- paste0("./models/MODEL_I_IRFs_", year, ".Rdata")
print(load(model_name))
var_irfs_cum_merged_2018 <- var_irfs_cum_merged

# Generate data for plotting 5-day responses to a permanent shock
final_input_2018 <- transform_var_irf(var_irfs = var_irfs_cum_merged_2018, days = 5)
plot_db_2018 <- make_plot_data(final_input = final_input_2018, 
                          which_effects = which_effects,
                          which_day = 5,
                          column_group = c("state_legislators"))
plot_db_2018$Year <- "2018"

final_input_2021 <- transform_var_irf(var_irfs = var_irfs_cum_merged_2021, days = 5)
plot_db_2021 <- make_plot_data(final_input = final_input_2021, 
                               which_effects = which_effects,
                               which_day = 5,
                               column_group = c("state_legislators"))
plot_db_2021$Year <- "2021"

plot_db <- rbind(plot_db_2018, plot_db_2021)
plot_db$rowgroup <- recode(plot_db$rowgroup,  "Trump" = "President")

# Make level variable
plot_db$level <- ifelse(plot_db$rowgroup %in% c("President", "Members of\nCongress",
                                               "National Media"), 
                        "National\ngroup", "State\ngroup")

#  Define the order for the plot
plot_db$direction <- factor(plot_db$direction,
                            levels = c("Row groups' response",
                                       "State legislators' response"
                            ))
plot_db$year_dir <- paste0(plot_db$direction, plot_db$Year, sep = "")
plot_db$rowgroup <- recode(plot_db$rowgroup,
                          "National Media" = "National\nMedia",
                          "State Media" = "State\nMedia",
                          "State Partisans" = "State\nPartisans")

# Generate Figure 
plot_db$level <- factor(plot_db$level, levels = c("State\ngroup", "National\ngroup"))

plot_overall <- ggplot(plot_db, aes(x = direction,  ymin = lwr, 
                                    ymax = upr, col = direction, pattern = Year)) +
  geom_hline(yintercept = 0, color = "red") +
  geom_linerange(aes(x=year_dir, ymin=lwr, ymax=upr), size = 5.5) + 
  geom_segment(aes(x=year_dir, xend=year_dir, y=pe-0.01, yend=pe+0.01), color="black", size=5) +
  geom_text(aes(x=year_dir, 
                y = upr+.23,
                label = Year), col = "black", size = 3.5) + 
  facet_nested(level + rowgroup ~ colgroup, switch = "y") +
  coord_flip() +
  ylab("\n5-day effect of a permanent 10 percentage point increase") +
  scale_color_manual("Effect direction:", values = c("gray80", "gray40")) +
  scale_shape_manual("Year:", values = c(2,1)) +
  guides(col = guide_legend(ncol=2,nrow=1, reverse = TRUE)) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 12),
    legend.title = element_text(size = 12),
    panel.background = element_blank(),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    strip.text.y.left = element_text(angle = 0, size = 11, 
                                     margin = margin(0,.2,0,.2, "cm")),
    strip.text.x = element_text(size = 11),
    axis.text.y = element_blank(),
    axis.title.x = element_text(size = 12),
    axis.text.x = element_text(size = 11),
    axis.ticks.y = element_blank(),
    axis.line = element_line(color = "black"),
    axis.title.y = element_blank()
  )
plot_overall

#  Save the figure
ggsave(paste0(root_dir, "/figures/figure02.png"), plot_overall, width = 8.5, height = 7)

